Poniżej przedstawiam kod MATLABa obliczający problem adwekcji-dyfuzji metodą elementów skończonych ze stabilizacją metodą minimalizacji reziduum, dla modelowego problemu Erikksona-Johnsona.
Przypomnijmy najpierw problem adwekcji-dyfuzji, problem modelowy Erikksona-Johnsona (opisany na przykład w dokumencie [1] zdefiniowany na obszarze kwadratowym \( \Omega = [0,1]^2 \) w następujący sposób: Szukamy funkcji \( u \) takiej że
\( a(u,v)=l(v) \forall v \) gdzie
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \)jest rozszerzeniem warunku brzegowego Dirichleta na cały obszar, natomiast
\( \beta = (1,0) \) reprezentuje wiatr wiejący z lewej strony na prawą, natomiast \( \epsilon = 10^{-2} \) oznacza współczynnik dyfuzji. Wielkość \( Pe=1/ \epsilon = 100 \) nazywa sie liczbą Pekleta, i definiuje ona wrażliwość problemu adwekcji-dyfuzji.
Kod można uruchomić w darmowym środowisku Octave.

Pobierz kod lub zob. Załącznik 5A.

W linii 830 podawany jest współczynnik dyfuzji \( \epsilon = 10 ^{-2} \). Jego odwrotnośc to tak zwana wartość Pecleta
\( Pe=\frac{1}{\epsilon} \).
Następnie tworzony jest wektor węzłów (zawsze liczby naturalne liczone od zera) \( knot\_trial\_x = [0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2] \) dla przestrzeni aproksymacyjnej,
oraz wektor odpowiadających im punktów wzdłuź osi x, z przedziału od 0 do 1, \( points\_x = [0 \quad 0.5 \quad 1] \) w których umieszczana będzie baza funkcji B-spline służaca do aproksymacji rozwiązania problemu Erikksona-Johnsona. Podobnie definiujemy wektor węzłów dla przestrzeni testującej \( knot\_test\_x = [0 \quad 0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2 \quad 2] \). Zgodnie z ideą metody minimalizacji reziduum, funkcji testujących musi być więcej niż funkcji aproksymujących. W naszym przykładzie używamy funkcji B-spline stopnia kwadratowego o ciągłości C1 dla przestrzeni aproksymacyjnej (trial) oraz funkcji B-spline stopnia trzeciego o ciągłości C1 dla przestrzeni testującej (test).
Wektory węzłów dla przestrzeni aproksymującej oraz testującej, i punkty na których rozpinamy wektory węzłów, definiujemy osobno dla osi x oraz dla osi y.
Kod może zostać uruchomiony w darmowym środowisku Octave.
Kod uuchamia się otwierając go w Octave oraz wpisując komendę
\( advection\_igrm\_adapt \)
Po chwili obliczeń kod otwiera dodatkowe okienko i rysuje w nim rozwiązanie numeryczne oraz dokładne. W drugim oknie kod rysuje reziduum, obliczane dodatkowo przez metodę minimalizacji reziduum, Jest ono miarą lokalnego błędu rozwiązania.
Kod oblicza błąd reziduum
\( Residual norm: L2 0.015808, H1 0.100417 \)
iKod oblicza również normę L2 i H1 z rozwiązania
\( Error: L2 45.90 procent, H1 61.85 procent \)
i porównuje do normy L2 i H1 z rozwiazania dokładnego.
\( Best possible: L2 2.02 procent, H1 14.20 procent \)
Widać że w celu poprawienia dokładności rozwiązania, konieczne jest zwiększenie siatki.

Autorzy kodów w MATLABie: Marcin Łoś i Maciej Woźniak.